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Introduction 


This paper studies a compact finite difference scheme for the solution of 
the elliptic system 


( 1 . 1 ) 


2 

div v - q u = 0 


v - P grad u = 0 


when u or its normal derivative are prescribed on the boundary of a domain. 
While important in its own right, this system also provides a prototype for 
the equilibrium equations for elastic materials. We shall show elsewhere that 
the methods described in this paper can be applied without essential 
modification to such equilibria problems. 

The compact finite difference scheme to be described represents a finite- 
volume method which expresses algebraic relationships between average values 
of u on the sides of a computational cell and the average values of the flux 
normal to the sides. The term "compact" refers to the fact that these 

relationships hold without reference to neighboring cells. An advantage of 
such schemes is that any prescribed boundary data may be incorporated with the 
same accuracy as the scheme itself if irregularly shaped cells are employed. 

Unfortunately, the algebraic problem presented by such schemes is 
difficult to treat, especially if fast iterative methods are sought. By a 
process analogous to eliminating the flux v_ in (1.1) it is possible to 
obtain algebraic relationships solely between the solution variables u in 
neighboring cells. This provides a two-stage iterative process, the first 
concerning cell neighbors in one direction, the second the neighbors in the 
other direction. These equations may be conveniently solved by a Gauss-Seidel 
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type of iterative process or, as will be shown elsewhere, by a multigrid 
method. Both stages may also be combined into one stage. 

We illustrate these features by treating first Laplace's equation 

2 

V u = 0. The combined one-stage method is seen to yield a second-order 
accurate nine-point finite difference scheme for the Laplacian. 

We conclude our discussion by describing the scheme for the more general 
problem (1.1) and indicating the sense in which the scheme is dissipative by 
means of an energy estimate. 

The method described in this paper has its origin in an earlier approach 
(Rose [1]) to which the reader is referred. 


2. A Compact Scheme for V u =■ 0 

Let (x ,y ) = (iAx.jAy), h = Ax/2, h = Ay/2 and let ir denote the 
i j x y x,j 

rectangular cell { |x-x | < h , ly-y.! < h }. We describe variables 

j y 

associated with the sides of tt , by referencing the center point of the 

i» J 

side; thus u^ ^ u^ j^ly^ indicate the average values of u associated 

with the sides of if. . 

i » j 

The translation operators s and t are defined by 


su i,j 5 "i+V 2 .r tu i,i s “i.j+v 2 

and we define central average and difference operators by 

y = (S + s -1 )/2 , y - (t + t _1 )/2 

x y 

6 = (s - s _1 )/2h , <$ = (t - t _1 )/2h . 

x x y y 
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Consider 
system form: 


the Dirichlet problem for 
if v - (v,w)', then 


V 2 u 


0 in a square domain. 


In 


( 2 . 1 ) 


div v = 0 


_v - grad u = 0. 


Corresponding to this system we consider the following compact finite 
difference scheme: in any cell if, u, v, w are related by the algebraic 

equations 

a) <5 v + <5 w=0 

x y 


( 2 . 2 ) 


b) p v - 6 u = 0, 
x x 


p w - 6 u = 0 

y y 


c) 


h^ 6 w 

y y 


(P - P )u = o. 
y x 


Equations (2.2a,b) are clearly consistent with (2.1). Equation (2.2c) 
2 

expresses an 0(h ) approximation to the value of u at the center of the 
cell; motivation for this particular approximation will be given in Section 6 
where energy estimates are discussed. We may expect this scheme to yield u 
to second-order and, noting (2.2b), (v,w) to first-order accuracy. 

Write 

U 5 ( W x u ’ 6 x u> M y U ’ 6 y ’ 


(2.3) 


X 

V = (p v, 8 v, P w 
x x y 
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Then (2.2) may be written in matrix form as 


(2.4) 


P V = Q U 


where 


( 0 1 0 

1 0 0 

0.0 1 

0 h 2 /2 0 

x 




0 0 0 

1 0 0 

0 0 1 

0-10 


If 


(2.5) 


R = P -1 Q 


the fluxes V may be expressed in terms of U by 


V = R U, 


( 2 . 6 ) 
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where, 


(2.7) 


" 0 

1 

0 

0 “ 


r r, “) 

-2 

0 

-2 

0 


— 1 

h 

-h 


—2 

0 

-2 

0 

0 

-2 

1 


-3 

-h 

0 

h 

0 


—4 

L 


and 


h 2 == (h 2 + h 2 )/2. 

V x yj 


. The value 


Now consider the contiguous cells it^ y ^ and v ^y^ j 

which is associated with the side common to these cells is given by 


(2.8) 


i.e • , 


V i,j + h x 5 x^ V i~ V2 . j 


= (y x " h x 6 xh+v 2 ,r 


(2.9) 


( (p - h 6 )s - (P + h 6 )s 1 )v. , » 0. 

' ' V V v' V V tr ' 11 


x xx x xx 


Recalling the definitions of P , <5 in terms of s as well as the 

X X 

definition of V in (2.3) we may write (2.9) as 


2 h x 0 x * - V 0, o)v = 0 
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so that, employing (2.6) and dropping the factor 2h x , we obtain, in terms of 
the rows Xl» Z2 


( 2 . 10 ) 


(S * it - “x ■ °> 


i »6 • j 


( 2 . 10 )' 


(h 2 6 2 - u (y - y )) u . 0. 
x xx y ' 


Similarly, by eliminating the value Wj common to the cells x i. 

/2 

and x we obtain 

i»j+ /2 


( 2 . 11 ) 


(h 2 S 2 - y (y - y )) u = 0. 
y y y x J 


For simplicity let us now assume h = h x = hy. Since 


.2 ,2 2 , 

h 6 = y - 1 , 

x x 


,2 .2 2 . 

h 6 = y - 1 

y y 


(2.10) and (2.11) assume the simpler form 


( 2 . 12 ) 


u. = y y u. , 
i.j x y i,j 


where u. * is a value on the side common to either x i , , , x , , 

1 »J V2 ,j i+ x /2 » j 

x ± j_ i y , x^ j + i^. Referring to the stencil indicated in Figure 1 


or 
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Figure 1 


(2.12) indicates that 


(2.13) 4u(P) = ulQj) + u(Q 2 ) + u(Q 3 ) + u(Q^,). 

Utilizing (2.12) to express u(Q) in terras of its neighbors as well we obtain 


(2.14) 


2 2 

u. ■ y v u. 
i» j x y i, j 


2 

which is a familiar second-order accurate nine-point expression for V u = 0. 

In the more general case to be described in a later section (2.13) appears 
in the more complex form 


(2.15) 


4 

u(P) = l u(Q i ) + u(P 1 ) + e 2 u(P 2 ). 


In this case the simplified discussion of iterative solutions of (2.12) to be 
given in the next section will, of course, not apply. Nevertheless the 
analysis is indicative of the more general situation. 
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3. Iterative Methods 

Here the structure of the difference equations obtained in the previous 
section is investigated in order to devise appropriate iterative methods for 
their numerical solution. 

We first rewrite (2.12) in the form 


(3.1) 


a) (1 - u x M y )» 1%j - 0, 

b) (* ' \ ’ °- 


Equations (3.1a) and (3.1b) represent difference equations for values of u 
defined respectively at the vertical and horizontal mid-points of the sides of 
the computational cells if. . . The problem is certainly fully determined 

1 9 J 

since at each grid point we have either an equation or a boundary value. 

We define a partition of the system (3.1) as follows: 

(i) Let p denote the vector of those unknowns defined at mid-points of 
vertical mesh lines. Equations (3.1a) represent finite difference 
equations defined at these points. 

(il) Let q denote the vector of these unknowns defined at mid-points of 
horizontal mesh lines. Equations (3.1b) represent finite difference 
equations defined at these points. 

With this partition (3.1) may be written In the following matrix form 



(3.2) 


9 
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where B Is an N(N-l) * N(N-l) matrix and I is the N(N-l) * N(N-1) 

identity matrix. The elements of B are negative and they contain at most 

four non-zero entries per row, these all being -1. The vectors b and c 
contain boundary values. 

The coefficient matrix of the system (3.2) is symmetric, positive definite 
and consistently ordered. The first two properties are important since then 
the convergence of the S.O.R. method is obtained for any value of the 
relaxation factor w in the range 0 < m < 2 (see Young [3]). The property 

of consistent ordering is important in the theory of the S.O.R. method because 

at present the calculation of the optimum relaxation factor is possible only 
for consistently ordered matrices. Results related to the determination of 
the optimum relaxation factor for the S.O.R. method are proved by Young [3]. 

Equation (3.2) illustrates the fact that each equation defined at a 
horizontal mid-point is not coupled to any other unknowns at horizontal mid- 
points. The same is true for equations defined at vertical mid-points. 

We note that this system is not equivalent to the one obtained by 
discretizing Laplace's equation using the rotated five-point formula on a 
uniform mesh. One reason for this is that the coefficient matrix of this 
system is reducible whereas that of (3.2) is not. 

Equation (3.1) may be written in component form as follows 

4p + Bq = b, 

(3.3) 

T 

B p + 4q = c. 

These equations have thus been cast in a form which is amenable to many of 
the standard iterative methods of solution. Equation (3.2) may be solved 
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using the following two-stage iterative process: 


4p 


(n+i) + Bq (n) = 


( (n+1) tn; 'i 

l u 4 » . ■ H V u', .J , 

^ i »j x y i 


(n) 


(3.4) 


B T p<" +1 > + 4, (n+1) » 


(„<»«> - U V u'"t 
^ l»j x y i, j 


(n+1) j 


This is the Gauss-Seidel method for solving the system (3.2). 

Each of the equations in the system (3.2) defines a rotated Laplacian 
operator acting on the value of u at the mid-points of the cells. The 
authors have had some experience with this operator while working with the 
Cauchy-Riemann equations. In that application serious problems with con- 
vergence were encountered. These difficulties appear to be due to treating 
only the first stage of a two-stage process. A similar effect is seen in 
A.D.I. methods for parabolic equations where the use of only one of the two 
steps leads to instabilities. 

As an alternative to the scheme (3.4) we may eliminate p from (3.3) to 
obtain the following nine-point scheme for q 


(3.5) 


( 161 - B T B)q = 4c - B T b 
i.e., Aq = d 


(u 




At interior points we have the following finite difference equation 

(3.6) 2 (q 1+1>J+1/2+ q 1>j+3/2 + 

+ ( q i+l,j+3/2 + q i+l,j- V 2 + ’i-l» j+3/2 + Vi, j- V 2 • ' 12 q i,J+V 2 ” °’ 
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at the point (x^»yj+ ly^ ) • We note that this is not the higher order nine- 
point approximation to Laplace's equation; rather (3.6) is second-order 
accurate. Boundary conditions are to be incorporated into interior equations 
in the usual fashion. The matrix A appearing in (3.5) is irreducibly 
diagonally dominant and so both the point Jacobi and point Gauss-Seidel 
iteration matrices are convergent and the associated iterative methods are 

convergent for any initial approximation q^ (see Varga [2]). This is 

evident, also, from the fact that (3.6) expresses q^ 1 j as a weighted 

average of its eight neighbors. 

It is interesting to note that although the coefficient matrix of the two- 
stage scheme (3.4) is consistently ordered, the coefficient matrix for the 

nine-point scheme (3.6) is not. An estimate of the optimum relaxation factor 
for the higher-order nine-point formula has been obtained by van de Vooren and 
Vliegenthart [4] using separation of variables. 

The rate of convergence of the two-stage scheme (3.4) is 0(h). 
Convergence may be accelerated by employing the S.O.R. method with optimum 
relaxation factor, in which case the rate of convergence increases to be 

0(h). 

When h * h or in the more general treatment of (1.1) to be described 
x y 

below the observation that the equations defined at the two sets of points are 
decoupled will no longer be true. In this case the equations corresponding to 
(2.8a) and (2.8b) will assume the more general form (2.15). This structure 
suggests the use of line relaxation as a method of solution; this would 
involve relaxing (2.8a) along horizontal lines and (2.8b) along vertical 


lines . 
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4. Numerical Results 

Here we demonstrate the performance of the two-stage method on a simple 
model problem. We consider Laplace's equation defined in the unit square with 
Dirichlet boundary conditions such that the true solution is given by 

u(x,y) = e ^ sin(iTx) . 

The S.O.R. method was applied to this problem with the optimum choice of 
relaxation factor. This choice of relaxation factor is given by 


1 + / [1 - a 2 ] 

where 

o = cost, — J . 

/2 

The reader is referred to Young [3] for the derivation of this formula. 

The algorithm was terminated when the magnitude of the maximum difference 
between successive iterates was less than 10~^. Table I shows the dependence 
of the relaxation factor, w, and the number of iterations on the mesh length. 

The behavior observed is typical of that of S.O.R. methods. We now have 
an effective relaxation method for solving the compact scheme in a form which 
should permit the development of a multigrid algorithm and thereby obtain an 
efficient means for solving the class of problems in this paper. We will 
report on this development separately. 
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Table I 


h 

(0 

No. of Iterations 

1/8 

1.57 

19 

1/16 

1.76 

34 

1/32 

1.87 

67 

1/64 

1.93 

129 


5. The General Problem 

2 

The discussion of V u = 0 given in Section 2 can be easily adapted to 

the more general system (1.1). 

Assume h = h„ = h and introduce the abbreviations 

A y 


(5.1) 


a = y 

x x 


P. 



P. 


Consider the compact scheme 


a) 

6 v + 6 w 

1 2 f 
-7 q 

y + 

y ] 

lu = 

0 

x y 

2 M v 

X 

y J 



b) 


y v 

- a 

6 

u = 

0 


X 

X 

X 



c) 


y w 

- a 

6 

u = 

0 


y 

y 

y 



d) 

h 2 (S w - I q 2 y u) - 
y 2 y J 

(a y 

- a 

y , 

lu = 

0 

v y y 

X 

X J 




As before, (a), (b), and (c) are clearly consistent with (1.1) and (d) is a 

n 

consistent center-point approximation of u to terms 0(h). 

Using the definitions of U, V given by (2.3) we may again express (5.2) 
in matrix form and, as a result, express V in terms of U by V = R U 
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where R Is the matrix described below. The result of eliminating v. . 

*■» J 

and w^ j at points common to neighboring cells is again of the form 


(5.3) 


a) (6 x £l - = 0 


b ) ( 6 y L 3 - % £ 4 ) U = 0, 


(c.f. (2.20)) in which _r^ is a row of R. 
From (5.2) we find 


(5.4) R = 


0 a 

(j q 2 + 0 

0 0 

-h" 2 o 0 

x 


-h~ 2 a 


a q 2 + h - 2 a J 


0 

0 

a 


Thus (5.3) leads explicitly to 


r , 2 x . 1,222 'i _ 

[ho o 6 - y a y - -r- h qy +y o y Ju = 0 

xxx xxx2 x xyy 


r 2 , p 1,222 'i n 

1 h o o 6 -y o y - — h q y +y o y u = 0, 
yyy yyy2 y yxx' 


or, using the identity 


2 

h (<5 0 6)u = (y 0 y)u - (y 0)u, 


and the definitions of o ,0 , to 

x’ y 
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(5.5) 


a) (y 2 p)u +~ q 2 h 2 u = y^ (P P)P y )u 

b) (y 2 p)u + 4 q 2 h 2 y 2 u = y ( (y p)y )u. 

2 y y^ x x' 


Reference to Figure 1 shows that eq. (5.5) involves 
center points of the sides of u. , i, . and , , i. . 

i±V2,j i>j±V2 

reduces to (2.10) and (2.11) when P = 1 and q = 0 

essential properties required to adapt to these equations 

methods described earlier. 


u at all of the 
Clearly, (5.5) 
and retains the 
the same iterative 


6 . 


An Energy Estimate 

Let 


Lu = - div p grad u + q 


2 


u 


and 



2 2 

J (p grad u + q u)dxdy. 
D 


Green's theorem applied to (1.1) yields, after multiplication by u, the 
familiar estimate for the energy norm Hull given by 


(6.1) Hull 2 = / u Lu dxdy + / uu ds. 

D r n 

In addition to providing a uniqueness argument for the solution of Lu = 0, 
(6.1) forms the basis for many important properties of the solution of this 
elliptic problem. 
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In this section we shall show that the compact difference scheme (5.2) (or 

(2.2) ) yields an identity closely related to (6.1). 

With respect to the operators y, S in a cell we shall write, if 
v = (v,w), 

(6.2) div, v = 6 v + 6 w, 

h — x y 

and also write the discrete form of Gauss' theorem as 

(6.3) \ div. v Ait = \ v*ji As. 

d h r 

Summation-by-parts is accomplished by the use of the identity 

(6.4) 6(<M>) = (y<f>)<Si|> + (y^p )64> . 


Finally, recall the definitions 



a = 

y 


y 


y 


p 


from (5.1) and set 

(6.5) L h u = - div h v + -| q 2 (u x + P y )«» 

(6.6) Hull 2 = l [o (6 u) 2 + o (6 u) 2 +j q 2 ( (u u) 2 + (y u) 2 )]Atu 

n d x x y y ^ x y 

We proceed as follows: first, multiply u by P x u and employ (5. 2d) 


to obtain 
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P u L u = - q 1 2 ( (u u ) 2 + (p u) 2 ) + h 2 (<$ w - t q 2 P u)‘ 

x ii 2^^x y y 2 n y J 


- (p u 5 v + P u S w). 
v x x y y 


Next, using (5.2b,c) (6.3), and (6.7), we may sum over cells in D to obtain 


(6.7) 


11 uH? + h 2 l (o /a )(6 w - 4 P u)^ AlT 


12 ^ 2 


x y' ' y 


) uu As + ), P u L.uAtt 

r n D x J 


so that 

(6.8) Hull 2 < y uu As + y p u L,uAtt, 

h r n D x " 

where the inequality is strict unless u = 0. 

This inequality implies the uniqueness and existence of the solution to 
(5.2). We leave it to the reader to adopt standard arguments to (6.8) to 

verify that the solution of the compact scheme provides a second-order 
approximation to (1.1) 

In order to help understand the effect of the dissipative term 

(6.10) h 2 y (o /a )(6 w - 4 q 2 P u) 2 Att 

P x y' v y 2 H y J 

in (6.7) again consider (2.2) where P = 1, q = 0. In place of (2.2c) 
consider 

1 2 

a h S w - (M - M )u = 0 

2 y y x 


(2.2c) 
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together with (2.2a,b). Let a denote the error between the solution of 
the modified compact scheme and the solution of the example discussed in 
Section 4. 

Table II compares II a # as a function of h and also of a. For this 
example the dissipative term (6.10) has the form 

a h 2 l (6 w) 2 Att, 

D y 

i.e., the dissipation is directly proportional to a. 



Values of He. II 
h,a 

Table II 

for different values of 

h and a 

a 

h - 1/8 

h = 1/16 

h = 1/32 

8 

0.0581 

0.0169 

0.0044 

1 

0.0123 

0.0030 

0.0009 

1/8 

0.0048 

0.0012 

0.0003 

1/16 

0.0042 

0.0012 

0.0009 

1/32 

0.0031 

0.0024 

0.0025 

0 

0.2798 

0.2858 

0.2871 


The effect of the dissipative term is clearly evident in the results given 

in Table II. For values of ot of 0(1) the convergence of the scheme is 
o 

0(h). However, for small values of a for which the dissipative terra 

2 

becomes less than 0(h) the convergence deteriorates. A particularly 

interesting feature occurs when a = 0, in which case the scheme fails to 
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converge. These results Indicate that the dissipative term Is required for 

o 

convergence of the scheme and should be 0(h ) in magnitude in order to 
2 

obtain 0(h) convergence. 

A closer examination of the transmission matrix for this example shows 
that R becomes singular for a + 0. As a result the compact scheme and the 
flux-elimination scheme are no longer equivalent in this limit. 


Conclusions 

We have described a compact system of finite difference equations for 
treating (1.1) and have shown how a related noncompact finite difference 
system provides an equivalent formulation. An energy estimate explains that 
the compact scheme is dissipative and also can serve to show that the solution 
approximates the solution of (1.1) to second-order accuracy. 

For V“u = 0 standard relaxation methods can be adopted to solve the non- 
compact scheme either as a two-stage method or, more directly, as a one-stage 
method. Both appear to be adaptable to multigrid solution methods. 

The methods described in this paper also apply, with little modification, 
to the equilibrium equations of elastic materials and appear to offer an 
interesting approach to convective-diffusion equations as well. We plan to 
report on these applications elsewhere. 
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